Derivation of transient relativistic fluid dynamics from the Boltzmann equation 
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In this work we present a general derivation of relativistic fluid dynamics from the Boltzmann 
equation using the method of moments. The main difference between our approach and the tradi- 
tional 14-moment approximation is that we will not close the fluid-dynamical equations of motion 
by truncating the expansion of the distribution function. Instead, we keep all terms in the moment 
expansion. The reduction of the degrees of freedom is done by identifying the microscopic time 
scales of the Boltzmann equation and considering only the slowest ones. In addition, the equations 
of motion for the dissipative quantities are truncated according to a systematic power-counting 
scheme in Knudsen and inverse Reynolds number. We conclude that the equations of motion can 
be closed in terms of only 14 dynamical variables, as long as we only keep terms of second order in 
Knudsen and/or inverse Reynolds number. We show that, even though the equations of motion are 
closed in terms of these 14 fields, the transport coefficients carry information about all the moments 
of the distribution function. In this way, we can show that the particle-diffusion and shear-viscosity 
coefficients agree with the values given by the Chapman- Enskog expansion. 

PACS numbers: 



I. INTRODUCTION 



Relativistic fluid dynamics is an effective theory to describe the long-distance, long-time dynamics of macroscopic 
systems, with important applications in relativistic heavy-ion collisions and astrophysics [l|. Relativistic fluid dynam- 
ics describes the conservation of (net) particle number and energy-momentum, 

d^N" = , dyT^ v = . (1) 

In general, these five equations contain 14 unknown fields, the four components of the particle 4-current iV M and the 
ten components of the (symmetric) energy-momentum tensor T^ v . Thus, these equations are not closed and one needs 
to specify nine additional equations of motion to solve them. The coefficients in the equations of motion (equation 
of state, transport coefficients, etc.) must be determined by matching fluid dynamics to the underlying microscopic 
theory. In the case of dilute gases, this is the Boltzmann equation. 

There are two widespread methods to provide additional equations of motion from the Boltzmann equation: the 
Chapman-Enskog expansion and the method of moments. In the Chapman-Enskog expansion [2|, the corrections 
to the single-particle distribution function in local equilibrium are assumed to be functions of the five traditional 
fluid-dynamical variables, temperature, chemical potential, and the three components of the fluid-velocity field, as 
well as gradients thereof. The corrections are systematically arranged in terms of an expansion in powers of the 
Knudsen number, given by the ratio of the mean-free path of the particles and a characteristic macroscopic length 
scale. As is well-known, the first-order truncation of the expansion leads to Navier-Stokes theory. Keeping second 
and higher-order terms one obtains the Burnett and super-Burnett equations, respectively Q. However, it has been 
shown that the Burnett equations suffer from the so-called Bobylev instability [J]. In the relativistic case, even the 
first-order equations, i.e., the relativistic generalization of the Navier-Stokes equations, are unstable @. 

Therefore, the relativistic extension of Chapman-Enskog theory should not be applied to derive the equations of 
relativistic fluid dynamics from kinetic theory. On the other hand, the method of moments [(| avoids the above 
mentioned problems. The method of moments was first developed by Grad [7[ for non-relativistic systems. In 
Grad's original work, the single-particle distribution function is expanded around its local equilibrium value in terms 
of a complete set of Hermite polynomials Q. This expansion is truncated and the distribution function is finally 
expressed in terms of 13 fluid-dynamical variables: the velocity field, the temperature, the chemical potential, the 
heat-conduction current, and the shear-stress tensor. In this case the heat-conduction current and shear-stress tensor 
become independent dynamical variables which satisfy partial differential equations that describe their relaxation 
towards their respective Navier-Stokes values. Grad's method is usually considered to be independent of the Chapman- 
Enskog expansion. However, we emphasize that Burnett-type equations can be obtained as the solution of Grad's 
equations in the long-time limit 0, [l0| ■ 
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Nevertheless, Grad's method has one major drawback: unlike the Chapman-Enskog expansion it lacks a small 
parameter, such as the Knudsen number, in which one can do power-counting and thus systematically improve the 
approximati on [111 ] . This deficiency, together with the bad performance of Grad's method in comparison to microscopic 
calculations [12|], have led researchers to abandon this approach for some time. However, recently a lot of effort has 
been made to reformulate the method of moments into a more reliable tool to describe nonequilibrium phenomena for 
large Knudsen numbers (l2| . For instance, in Ref. fl3l | Grad's equations were regularized to have a wider domain of 
validity in Knudsen number and then shown to be in good agreement with microscopic calculations. Such approaches, 
however, were only formulated for non-relativistic systems. 

The generalization of Grad's method of moments to relativistic systems has been pursued by several authors [T3 |. 
The most widely employed approach is due to Israel and Stewart [Ifj. Here, the distribution function is expanded 
around its local equilibrium value in terms of a series of (reducible) Lorentz tensors formed of particle 4-momentum 
k^, i.e., 1, k^, k^k u , .... In Israel and Stewart's 14-moment approximation one truncates the expansion at second 
order in momentum, i.e., one only keeps the tensors 1, k^ , and k^k v , with 14 unknown coefficients (the trace of k^k 11 
is equal to m 2 , the rest mass of the particles) to describe the distribution function. The coefficients of the truncated 
expansion can then be uniquely related to the 14 components of the particle 4-current, N^, and the energy-momentum 
tensor, T^ v , the so-called matching procedure. While particle and energy-momentum conservation ([1} are obtained 
from the zeroth and the first moment of the Boltzmann equation, the additional nine equations of motion follow from 
the second moment of the Boltzmann equation. However, Israel and Stewart's theory shares the same problems of 
Grad's original approach: it lacks a parameter in which one can do systematic power-counting of corrections to the 
local equilibrium distribution function. 

It was recently confirmed that, at least for some special problems, the Israel-Stewart equations [l5| are not in good 
agreement with the numerical solution of the Boltzmann equation fl6l fl7j . Initial attempts to improve Israel and 
Stewart's theory were already made in Refs. [l8l - [2Cl |. but Israel and Stewart's 14-moment approximation was still 
used. In this paper we demonstrate that Israel-Stewart theory, as well as all previous attempts to improve it are 
actually incomplete. The reason is that the 14-moment approximation neglects infinitely many terms of first order 
in the Knudsen number. In our approach all terms of the moment expansion are included and the exact equations 
of motion for these moments are derived. These exact equations still contain the degrees of freedom and microscopic 
time scales of the Boltzmann equation. We prove that, in order to derive a causal dynamical equation for a given 
dissipative current, it is necessary to resolve at least the slowest corresponding microscopic time scale arising from 
the Boltzmann equation, in agreement with the results of Ref. [Toj . Unlike in Israel- Stewart theory, the truncation 
of the resulting equations of motion in terms of only 14 dynamical variables is then implemented by a systematic 
power-counting scheme in Knudsen number, Kn, and in the ratios, Rjj 1 = |IT| /P , R^ 1 = |n M | /n , R^ 1 = /Po, 
where n is the bulk viscous pressure, n M is the particle-diffusion current, n^ 1 ' is the shear-stress tensor, and Po and 
no are the pressure and the particle density in local equilibrium, respectively. The ratio R~* is related to the inverse 
Reynolds number in non-relativistic situations. We shall in somewhat loose terminology refer to all of them as "inverse 
Reynolds numbers" in the following. The resulting fluid-dynamical equations and coefficients are different from the 
ones obtained via the 14-moment approximation. We calculate the numerical values of the coefficients for a massless 
classical gas. We show that our values for the heat-conductivity and shear-viscosity coefficient agree with the ones 
calculated via Chapman-Enskog theory [6j. 

This paper is organized as follows. In Sec. [II] we review how fluid-dynamical variables are extracted from the Boltz- 
mann equation. In Sec. IHII we demonstrate how to expand the single-particle distribution function in terms of a com- 
plete, orthogonal basis in momentum space. In contrast to Israel and Stewart's non-orthogonal basis 1, fc M , k^k 1 ' , . . ., 
our approach uses irreducible tensors in 4-momentum fc M , and is thus orthogonal. The coefficients of the irreducible 
tensors in the expansion of the single-particle distribution function are orthogonal polynomials in the rest-frame en- 
ergy and moments of the correction to the equilibrium distribution function. Section IIVI derives an infinite set of 
equations for these moments, which is still completely equivalent to the Boltzmann equation. In Sec. [V] we introduce 
our power-counting scheme in terms of Knudsen and inverse Reynolds numbers. Then, by diagonalizing the linear part 
of the set of moment equations, we demonstrate how to identify the slowest microscopic time scale of the Boltzmann 
equation for each dissipative current. We shall derive dynamical equations for the slowest modes, but approximate 
faster modes by their asymptotic solution for long times. This will then lead, in Sec. I VII to the complete set of 
fluid-dynamical equations which contains all terms up to second order in Knudsen and inverse Reynolds numbers, i.e., 
C(Kn 2 , R^RJ" 1 , KnR" 1 ). In Sec. lVIII we first demonstrate the validity of our approach by restricting the calculation 
to the 14-moment approximation and recovering the results of Ref. [20( for the transport coefficients for the case of an 
ultrarclativistic, classical gas with constant cross section. We then show how to successively improve the expression 
for the transport coefficients by extending the number of moments to 14 + 9 x n. We explicitly study the cases n = 1, 2, 
and 3. We end this work with a discussion and conclusions in Sec. I Villi Various appendices contain intermediate 
steps of our calculations. We use natural units K = c = kg = 1- The metric tensor is = diag (+, — , — , — ). 
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II. FLUID-DYNAMICAL VARIABLES FROM THE BOLTZMANN EQUATION 

We start with the relativistic Boltzmann equation, 

fc^/k =C[f], (2) 



where — (fc°, k), with k° = Vk 2 + to 2 and m being the mass of the particles. For the collision term, we consider 
only elastic two-to-two collisions with incoming momenta fc, fc', and outgoing momenta p,p', 

C{f] = ^J dK'dPdP'W kk ^ pp , (/p/p>/k/k< - /k/k'/p/p') , (3) 

where v is a symmetry factor (= 2 for identical particles), Wkk/->pp/ is the Lorentz- invariant transition rate, and 
dK = gd 3 \i/ [(27r) 3 fc°] is the Lorentz-invariant momentum-space volume, with g being the number of internal degrees 

of freedom. We introduced the notation /k = 1 — afk, where a = 1 (a = —1) for fermions (bosons) and a = for a 
classical gas. 

In kinetic theory, the conserved particle current iV M and the energy-momentum tensor T" u are expressed as moments 
of the single-particle distribution function, 

= (k") , T» v = (kW) , (4) 

where we adopted the following notation: 



= fdK (■■■)&. (5) 



The particle current and the energy-momentum tensor can be tensor-decomposed with respect to the fluid 4- velocity 
w M . To this end, we have to specify the rest frame of the fluid. We introduce u" as a time-like, normalized (u M it M = 1) 
eigenvector of the energy-momentum tensor, 

T» v u v = eu", (6) 

where the eigenvalue e is the energy density, i.e., we work in the Landau frame pH ). Next, we divide the momentum 
of the particles k" into two parts: one parallel and one orthogonal to u M , 

fe M = E k u" + A>> , (7) 

where we defined the scalar E^ = u^k" = u ■ k and used the notation — A^A" , with A" u = g" v — u^u u being 
the projection operator onto the 3-space orthogonal to it M . 
Then, the tensor decomposition of N" and T' 1 " reads 

N" = nil" + ?A T" u = e u"u lJ - A"" (P a + II) + tt^ , (8) 

where the particle density n, the particle-diffusion current n M , the energy density e, the shear-stress tensor 7r M,y , and 
the sum of thermodynamic pressure, Pq, and bulk viscous pressure, II, are defined by 

n=(E k ) , n" = ^>) , e=(El) , tt"" = , P + H = -i (A^k^K) , (9) 

where A^") = A^A al3 and = [A£A£ + A^A^ - (2/3) A^A Q/3 ]/2 denotes a projector onto that part of a rank-2 
tensor, which is symmetric, orthogonal to u" , and traceless. 

Next, we introduce the local-equilibrium distribution function as /ok = [exp (/3q E^ — c*o) + a] , where /3q and «o 
are the inverse temperature and the ratio of the chemical potential to temperature, respectively. The values of ckq 
and fto are determined by the matching conditions, 

n = n = (£ k )o, £ = £ o = (-££)„ » (10) 

where 

(•••)o= / dK(-..)fou. (11) 



Then, the separation between thermodynamic pressure and bulk viscous pressure is achieved as Po — — (A^"fc Al fc^) /3 
and n = - {A^k^k v ) s /3, where 

(•••>* = <"->-<"->0. (12) 

The fluid-dynamical conservation laws ([1]) are equations of motion for n, e, and u M , hence one needs nine additional 
equations to determine the dissipative corrections II, n^, and tt^ 1 ". In the following, we shall use the method of 
moments to derive these equations. 
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III. EXPANSION OF THE SINGLE-PARTICLE DISTRIBUTION FUNCTION IN TERMS OF 

IRREDUCIBLE TENSORS 



In this section, we expand the single-particle distribution /k in terms of irreducible tensors. It is convenient to 
factorize the local-equilibrium distribution function /ok from /k, 



h = /ok (1 + ./ok0kj , (13) 

where 0k represents the deviation from local equilibrium and is a function of x^ and k^, which is ultimately determined 
by the solution of the Boltzmann equation @. 

The next step is to expand 0k in terms of a complete basis of tensors formed of k^ and E^. As mentioned in the 
introduction, Israel and Stewart chose the following basis to expand 0k : 1, k 11 , k^k u , fc''fc I/ fc A , . . ., and then truncated 
the expansion after the second-rank tensor k^k v , that is 0k = £k + e k^t + ^ k^kv, where ek, e^, are ^ ne expansion 
coefficients [l5j . Note that these tensors are not irreducible with respect to Lorentz transformations A^ v that leave 
the fluid 4-velocity invariant, A^ v u v — u^. As a consequence, they are also not orthogonal, see Chapter VI, Sec. 
2a of Ref. @. Therefore, the expansion coefficients cannot be straightforwardly obtained: in a non-orthogonal basis, 
this requires in general the inversion of an infinite-dimensional matrix. Also, this implies that the exact form of 
the expansion coefficients cannot be obtained once the expansion is truncated. Therefore, the approach of Israel and 
Stewart does not provide the complete expressions for the expansion coefficients. 

In order to avoid such problems, we expand 0k using the irreducible tensors, 

I, fc<"> , k^k v) , k^k v k x) (14) 

as a basis. It should be emphasized that these tensors form a complete and orthogonal set, analogous to the spherical 
harmonics [22]. These irreducible tensors are defined by using the symmetrized and, for m > I traceless, projection 
orthogonal to u M as 

A <m~Mm) = A%:::%A vi ""' m , (15) 

where the projectors A^ 1 .'.'.^™ are defined in Ref. @, see Appendix [Fl for details. In order to obtain the irreducible 
tensors (fT4)) , we apply the projection (IT5t to A Ul "' Vm = k Ul ■ ■ ■ k Vm . The tensors (fT4|) satisfy an orthogonality condition, 



J dKF k k^ ■ ■ ■ fc"»> k {ui ■ ■ ■ k Vn) = (2 ^ m 1 " ) ,, A^;;;C / dKF k (A^k a k ) m , (16) 



where n, m = 0,1,2,.. ., is an arbitrary function of and S mn denotes the Kronecker-delta. Using the basis (|14p . 
0k can be expanded as 



oc 

9k = 



J^xt-^k^-.k^. (17) 



The index I indicates the rank of the tensor A^ 1 and i = corresponds to the scalar A. The coefficients A k M 

(i) 

are complicated functions of Ek and are further expanded in terms of an orthogonal basis of functions P^' , 



A < ri -M = ^ c ( r ,.. w)p w ) (18) 



n=0 



where Ni is the number of functions P^n considered to describe the £-th rank tensor A^ 1 . In principle, Ni should 
be infinite, however in practice, the expansion (|I8p must be truncated and Nt characterizes the truncation order. The 



function P kn are chosen to be polynomials of order n in energy, Pk, 



M _ 

r=0 

which are constructed to satisfy the orthonormality condition 



^=E^k, (19) 



^- W « = U (20) 
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where oj^' is defined as 



— (A^k a k ) /ok/ok. (21) 



The coefficients an} and the normalization constants can be found via Gram-Schmidt orthogonalization using the 
orthonormality condition (|2U)) . see Appendix [E] for details. We note that, in the limit of massless, classical particles, 

the polynomials P^ n correspond to the associated Laguerre polynomials. 

Since the expansion (fl7|) employs an orthogonal basis, the expansion coefficients in Eq. (|T8|) can be immediately 
determined using Eqs. (fT5|) and (j2"0)) . For n < Ng they are given by 

For the sake of later convenience, these expansion coefficients are re-expressed as linear combinations of irreducible 
moments of <5/ k = /k - /ok, 

pJt 1 -"* = (e£ k^ 1 ■■■k t " ) ) s , (23) 



such that 



^^E^ 1 "*. ( 24 ) 

n=0 



where we defined the energy-dependent coefficients 



(£) = A (fl 



m—n 

Consequently, the distribution function itself can be expressed as a series in the irreducible moments (|23j) of <5/k, 

oo Nt 



/k = /ok + /ok/ok^^<^ 1 -^A :(ui - . (26) 



£=0 n=0 

We remark that the matching conditions and the definition of the velocity field imply that p\ = p2 = Pi — 0. 



IV. GENERAL EQUATIONS OF MOTION 

The time-evolution equations for the moments pjf 1 ' " Mf can be obtained directly from the Boltzmann equation by 
applying the comoving derivative to the definition (|23|). together with the symmetrized traceless projection, 

p<w-w> = A£;;;« J dKElk^ ■ ■ ■ k ^6f k , (27) 

where A = u^d^A = dA/dr and pf 1 ' = A{£.'.'.'# . Using the Boltzmann equation © in the form 

<5/k = -/ok - E^KVfow - E^k u V v 5f k + E^C [/] , (28) 

where V M = A^c^, and substituting this expression into Eq. (1271) . one can obtain the exact equations for the comoving 
derivatives of 

Using the power-counting scheme developed in Sec. |Vl we will show that, in order to derive the equations of motion 
for relativistic fluid dynamics, it is sufficient to know the time-evolution equations for the moments (|23p up to rank 
two, i.e., for p ri p%, and p% u . Similar equations could also be derived for higher-rank irreducible moments, if needed. 
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Thus, using Eqs. ([27]) and ([28]). we obtain 



C r _x = 4°) 6 - ^-W + ^tt""^ + ^d^n" + (r - 1) p£V^ + rp^ - 



Do 



Do 



-- [(r + 2)p. r -(r-l)m 2 p r _ 2 ] 0, 
p>> - = a«I* + p>t + i [(r - 1) m 2 p?_ 2 - (r + 3) - A^p*^ + rp^K 
+ ~ [(2r - 2) m 2 p*_ 2 - (2r + 3) p v r ] + 1 [m 2 rp r _x - (r + 3) p r+x ] 

+ _ yAln + _ i v „ ^ m 2 pr i _ pr+i) + (r _ i) 

£o + "o <J 
• _ C (M = 2a ( 2)(J „ _ 2 J (2r + 5) ^ _ m , 2 (j . _ 1} p A(M j ^) + 2 ^ w ,) 

+ ^ [(r + 4) p r+2 - (2r + 3) m 2 p r + (r - 1) m 4 pr-- 2 ] a*" + j|V^ (p^ - m a p r -i) 
2 



(29) 



(30) 



(r + 5) pf +1 rm 2 ^ u"> - - [(r + 4) pT - m 2 (r - 1) p^ 2 



+ (r - 1) p;_ 2 ^<7 X p - A^pVxP r -l + r Pr-l u * . 

where we introduced the generalized irreducible collision terms 

C^ 1 "^ = [ dKE^k^ 1 ■■■k^ t) C[j] 



(31) 



(32) 



We further defined the shear tensor a^ v = V^u v \ the expansion scalar 9 = V M it M , the vorticity tensor uj^ = 
{y y, u v — V"/) /2 and introduced 1^ — V^ao- All comoving derivatives of ocq and /3 that appeared during the 
derivation of the above equations were replaced using the exact equations obtained from the conservation laws of 
particle number, energy, and momentum, 





1 


&Q = 


D^ 




1 


$0 = 






1 


it* 1 = 




£o + 



The coefficients af°\ qJ- 1 ', and a>r ) are functions of temperature and chemical potential and have the general form, 



a = — {-J 3a {n 9 + d^) + J 20 [(e + Po + II) 9 - n^a^}} , 
V, = 7^ { --ho {n + d^) + J 10 [(e + P + H) 9 - n^a^]} , 
(V^Po - Ylu" + VII - A^tt q ^) . 

--10 

.,(2) 



(33) 
(34) 
(35) 



(0) 



(1 -r) J r i - Zro 



1 



D 



20 



[G 2r (e + Po) - G 3r n ] 



n 



+2,1, 



4 = Jr+1,1 

£0 + "0 

4 2) = J r+2,1 + (r - 1) 4+2,2 , 



where we defined the thermodynamic functions 
I m (a ,A))= ( 2 g |i)M 



3/„ 



A) 



Jn^JmO Jn— 1,0*^771+1,0 ) D n q — Jn-\-l.qJn—l,q J ri 



(36) 
(37) 
(38) 

(39) 
(40) 



The dissipative quantities appearing in the conservation laws can be (exactly) identified with the moments 



Po 



£.11, pS = n", p ^ = ^ 



(41) 



We note that the derivation of these general equations of motion is independent of the form of the expansion of the 
single-particle distribution we introduced in the previous section. 
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V. POWER COUNTING AND THE REDUCTION OF DYNAMICAL VARIABLES 



So far, we have derived a general expansion of the distribution function in terms of the irreducible moments of <5/k, 
as well as exact equations of motion for these moments. There is an infinite number of equations (labeled by the 
index r), and the equations for the moments up to rank two, Eqs. (|29p - (|31|) . contain moments of rank higher than 
two. In general, one would have to solve this infinite set of coupled equations in order to determine the time evolution 
of the system. However, in the fluid-dynamical limit, it is expected that the macroscopic dynamics of a given system 
simplifies, and therefore it can be described by the conserved currents and T^ v alone. 

From the kinetic point of view, it is usually assumed that the validity of the fluid-dynamical limit can be quantified 
by the Knudsen number, 



where £ m i CT and L macr are typical microscopic and macroscopic length or time scales of the system, respectively. 
The relevant macroscopic scales are usually estimated from the gradients of fluid-dynamical quantities, while the 
microscopic scales are of the order of the mean-free path or time between collisions. It is generally assumed that when 
there is a clear separation of the microscopic and macroscopic scales, i.e., when Kn <C 1, the microscopic details can 
be safely integrated out and the dynamics of the system can be described using only a few macroscopic fields. 

Furthermore, we also expect fluid dynamics to be valid near local thermal equilibrium, i.e., when 5fk <C fok- We can 
quantify the deviation from equilibrium in terms of the macroscopic variables by defining a set of ratios of dissipative 
quantities to the equilibrium pressure or density. These can be understood as generalizations of the inverse Reynolds 
number and will be denoted as 

R-.H.H-.M.B-.M. (43, 
■mj n o "o 

Since the non-equilibrium moments are integrals of Sfk while the equilibrium pressure and particle density are integrals 
over the equilibrium distribution function /ok, these ratios quantify the deviations from equilibrium. 

With this in mind, it is clear that these two measures, the Knudsen number and the inverse Reynolds number, 
can be used to quantify the proximity of the system to the fluid-dynamical limit. In general, these two measures 
are independent of each other, e.g. a system can be initialized in such way that the Knudsen number is large, but 
the inverse Reynolds number is small or vice versa. When deriving transient fluid dynamics, one should not a priori 
assume that Kn ~ Ri~ 1; while the Reynolds and Knudsen numbers are certainly related, their relation is in principle 
dynamical and is precisely what we aim to find. Only for asymptotically long times, the solutions of the dynamical 
equations yield Kn ~ R;" 1 ; as will be discussed in more detail below. 

In the traditional 14- moment approximation introduced by Israel and Stewart [l5j . the fluid-dynamical limit is 
implemented by a truncation of the expansion of the distribution function, which corresponds neither to a truncation 
in Knudsen nor in inverse Reynolds number. In this sense, the domain of validity of the equations of motion obtained 
via the traditional I4-moment approximation is not clear, because it is not possible to determine the order of the 
terms that were neglected. In order to obtain a closed set of macroscopic equations with a clear domain of validity 
in both Kn and R" 1 , another truncation procedure is necessary. The derivation of this is the main purpose of this 
section. 

First, we re-write the collision terms w ^ by linearizing the collision operator C[f] in the deviations from the 

equilibrium distribution functions. We then use the moment expansion (I26|) to obtain 
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micr 



(42) 



Kn = 



L 



macr 




(44) 



where 




x 



fjr-lu(vi ...uvi) f't/Wh, ■ ■■k v -I- "W W k', ■ ■ ■ k' ■ • • n « - n', ■ ■ ■ rl > 

K K ^kn^^i K vi) + rt k'n K (Ki K vt) n pnP(vi Put) n p>nP(vi P v t ) 



) ■ (45) 



The details of the derivation are relegated to Appendix [S] The coefficient A r7 l is the (rn) element of an (Nt + 1) x 
(Ni + 1) matrix and contains all the information of the underlying microscopic theory. We remark that, for £ = 0, 
the second and third rows and columns (r, n — 1,2) and, for i = 1, the second row and column (r, n = 1) are zero, 
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because the moments pi, P2, and p± vanish due to the definition of the velocity field and the matching conditions, 
Eqs. © and (|10l) . Therefore, in order to invert A^\ for I — 0, we have to exclude the second and third rows and 
columns and, for I = 1, the second row and column. 

As already mentioned, fluid dynamics is expected to emerge when the microscopic degrees of freedom are integrated 
out, and the system can be described solely by the conserved currents. The exact equations of motion ([29]) - (|3Tj) 
contain infinitely many degrees of freedom, given by the irreducible moments of the distribution function, and also 

infinitely many microscopic time scales, related to the coefficients Am- As was argued in Ref. [loj . the slowest 
microscopic time scale should dominate the dynamics at long times, i.e., in the fluid-dynamical limit. In order to 
extract the relevant relaxation scales, we have to determine the normal modes of Eqs. (1291 - (f5T|) . i.e., we diagonalize 
the part which is linear in the irreducible moments p^}'"^ 1 . These are the linear terms on the left-hand sides arising 
from Eq. (|44l) and the first terms on the right-hand sides. The nonlinear terms from Eq. (1441) as well as the remaining 
terms on the right-hand sides, which are nonlinear in the moments or are gradients of moments, are not considered in 
the diagonalization procedure. Identifying and separating the microscopic time scales of the Boltzmann equation is 
also the basic step for obtaining general relations between the irreducible moments and the dissipative currents and, 
as we shall see, closing the equations of motion in terms of iV M and T^ v ' . 

For this purpose, we shall introduce the matrix fiW which diagonalizes A^\ 

(0- 1 ) W AW0W=diag( X «,...,xf,...), (46) 



where xf are the eigenvalues of A^\ We further define the tensors X^ 1 '"^ as 

N, 

,-nC) 



j=o 

These are the eigenmodes of the linearized Boltzmann equation. Multiplying Eq. (l4"4l with (Q -1 )^ from the left and 
using Eqs. (|4"6")l and (|4"T)) we obtain 

C j-i " W> = -xf )x f "' W + ( terms nonlinear in Sf) . (48) 

3=0 

where we do not sum over the index i on the right-hand side of the equation. Then we multiply Eqs. (|29p - pip with 
(rr 1 )^ and sum over r. Using Eq. (|48|) . we obtain the equations of motion for the variables xj 11 '"^ 1 , 

Xi + xf ) X l = f3f ] 6 + (higher-order terms) , 
X^ ] + x { pX? = + (higher-order terms) , 

X^ v) + x^Xr = P- 2) <J^ + (higher-order terms) , (49) 

where we introduced the coefficients 

N Ni N 2 

/f 3 = E r 1 )^ ^= E (n- 1 )£ ) «S 1) . ^ 2) = 2 E(^ 1 )l 2)a f- ( 5 °) 

j=0,#l,2 J=0,#1 3=0 

With "higher-order terms" in Eqs. (|49p we refer to the terms nonlinear in Sf from Eq. (|48p as well as to the nonlinear 
and gradient terms on the right-hand sides of Eqs. (|29p - P"Tj) . As expected, the equations of motion for the tensors 
X^x—iii decouple j n the linear regime. Without loss of generality, we order the tensors X£ ll "' w according to increasing 

xi l \ e.g., in such a way that xi < Xr+n ^ 

By diagonalizing Eqs. ((2^)) - PT|) we were able to identify the microscopic time scales of the Boltzmann equation 

given by the inverse of the coefficients Xr ■ It is clear that, if the nonlinear terms in Eqs. (|49|) are small enough, each 
tensor X^ 1 '" 1 - 11 relaxes independently to its respective asymptotic value, given by the first term on the right-hand 
sides of Eqs. (|4"9")l (divided by the corresponding Xr^), on a time scale ~ 1/xr ■ We will refer to these asymptotic 
solutions as Navier-Stokes values. By neglecting all these relaxation scales, i.e., taking the limit Xr — * oo with 
fir^ /x^ fixed, all irreducible moments p^,' 1 ' " w become proportional to gradients of oto, /3q, and u M , and we obtain a 
Chapman-Enskog-type solution, which at first order in the Knudsen number results in the relativistic Navier-Stokes 
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equations of fluid dynamics. As already mentioned in the introduction, this type of solution is unstable and acausal, 
hence it cannot serve as a proper description of relativistic fluids. 

The solution for this problem was also mentioned in the introduction. To obtain causal and stable equations one 
must take into account the characteristic times within which the bulk viscous pressure, the particle-diffusion current, 
and the shear-stress tensor relax towards their asymptotic Navier-Stokes values. As shown in Ref. [lOj, in the fluid- 
dynamical limit these are given by the slowest microscopic time scales of the underlying microscopic theory, i.e., the 
fast relaxation scales are not expected to contribute. 

In practice, this is implemented by assuming that only the slowest modes with rank 2 and smaller, Xq, Xq, and 
Xq , remain in the transient regime and satisfy the partial differential equations (1491) . 

Xo + x °^ Xo = + (higher-order terms) , 
X^ ] + Xu ] Xu = Ppl" + (higher-order terms) , 
X^ V) + Xo ] Xg v = (3 { 2) a^ + (higher-order terms) , (51) 

while the modes described by faster relaxation scales, i.e., X r , X^, and X^ v , for any r larger than 0, will be 
approximated by their asymptotic solutions, 

X r ~ — 77y\® + (higher-order terms) , 

Xr 

X^ ~ —(T)!^ + (higher-order terms) , 

Xr 

B (2) 

X? v ~ -Wc^ + (higher-order terms) . (52) 

Xr 

While this approximation is similar to the Chapman-Enskog expansion, Eqs. (|5I[) go beyond the Chapman-Enskog 
expansion by including the transient dynamics. 

Note that, for r > I, X r , X^, and X^ v are of first order in Knudsen number, O(Kn). The reason is that the 

gradient terms 9, I^, and a^ v are proportional to L^^,., while 1/xr is proportional to 4nicr- The coefficients 
are simply functions of the thermodynamic variables ao, /3q> and thus of order 0(1). 

Furthermore, in order to obtain the traditional equations of fluid dynamics given in terms of the conserved currents, 
there should not appear any tensor X^ vX - - with rank higher than 2. Neglecting such tensors can be justified by proving 
that they have asymptotic solutions which are at least C(Kn 2 , KnR^ 1 ), i.e., beyond the order we consider here. 

Equations (|52[) enable us to approximate the irreducible moments that do not appear in the conserved currents in 
terms of those that do occur, namely the particle-diffusion current, the bulk viscous pressure, and the shear-stress 
tensor. We now show how to do this. We first invert Eq. (|4T|) . 

= , (53) 

3=0 



then, using Eqs. (|52|) . we obtain 



No o(0) 



Pi ~ «£> * + £ Og> % 9 = Sl$X + O(Kn) , 

j=3 Xj 
Ni 



Pi 



>» ~ n$XtT + J2 ^f%j^ = nftxp + 0(Kn) . (54) 



3 =2 Xj 

N 2 o(2) 
I P 3 



(2) 

3 = 1 X) 



Here, we indicated that the contribution from the modes X r ,X^, and X^ v for r > 1 is of order O(Kn). 

Taking % = in the above equations and, without loss of generality, setting Oqq = 1, we obtain from Eqs. (|4I|) the 
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relations 



3=3 A 3 

^0 - n Z^"0j (i) J > 

3=2 Xj 



^2 o(2) 

Xr^^-^og^a^. (55) 

3 = 1 Xj 



Substituting Eqs. ([35) into Eqs. ([55]). 



m 2 



3 



0(Kn 2 , KnR" 1 ) . (56) 



To obtain Eqs. (156]) , we further used that X^ 1 '"^ ~ C(Kn 2 ,KnR i *) for ^ > 3, and defined the transport coefficients 

Ci = ^ ] T ir a r ^ J K 7ii — ^ ] 7"i r a r ^ i % = ^ ] T ir a r ^ > 

where we introduced the inverse of ,4«, tW = (y^ 1 )^' and used the relation, 

r m - 2^ ™ (I) \ >mn ■ 

In the next subsection, we shall identify the coefficients Coi K n0, and rj as the bulk- viscosity, particle-diffusion, and 
shear-viscosity coefficients, respectively. 

So far we have proved that, by taking into account only the slowest relaxation time scales, all irreducible moments 
^a> Q £ ^ e deviation of the single-particle distribution function from the equilibrium one can be related, up to first 
order in Knudsen number, O(Kn), to the dissipative currents, II, n M , and 7r Miy . This demonstrates that in this limit, 
it is possible to reduce the number of dynamical variables in Eqs. (|29[) - (1311) to quantities appearing in the conserved 
currents. This will be explicitly shown in the next section. 

We remark that similar relations between the irreducible moments and the dissipative currents can also be obtained 
with the 14-moment approximation, but with a different set of proportionality coefficients. However, in the traditional 
14-moment approximation such relations are obtained by explicitly truncating the moment expansion (|26l) and, as 
a result, they are not of a definite order in powers of Knudsen number. This is the reason why the 14-moment 
approximation does not give rise to equations of motion with a definite domain of validity in Knudsen and inverse 
Reynolds numbers. 

Note, however, that the relations (151)1) are only valid for the moments p^ vX " with positive r. This is not a problem 
since similar relations can also be obtained for the irreducible moments with negative r. We expect the expansion 
(|26[) to be complete and, therefore, any moment that does not appear in this expansion must be linearly related to 
those that do appear. This means that, using the moment expansion, Eq. (|26[) . it is possible to express the moments 
with negative r in terms of the ones with positive r. Substituting Eq. (|26[) into Eq. (|23p and using Eq. (|16[) . we obtain 

p-v^E^Vr^ (58) 

n=0 

where we defined the following thermodynamic integral 

H$ = ^yyyy / dK fojo^n^ (A^k a k^ . (59) 
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Therefore, Eqs. ([561) lead to 



where we introduced the coefficients 



P-r 


= __L~(o)- 

2 IT 




m 


P-r 


= 7< 1 >n" + < 


LLV 

P-r 


= 4 2 K^ + 


\XV" ■ 

'— r 


= C(Kn 3 ) , 



0(Kn) , 



(60) 



No Ni N 2 



„(o) _ v ^°)n (0) -yW- V j-«fi (1) ^(2) - Vjc-( 2 )n (2) f6D 

/r ~~ / y •'rn "nO ' /r ~~ / , • r rn i!: nO ' fr — / t • r rn s L n0 • V u± ^ 



VI. COMPLETE FLUID-DYNAMICAL EQUATIONS TO SECOND ORDER 

Now we are ready to close Eqs. (J2S|) - ([3"Tj) in terms of the dissipative currents appearing in and and derive 
the fluid-dynamical equations of motion. For this purpose, it is convenient to use the inverse 
which naturally satisfies t^'A^' = 1. Hence, it is straightforward to rewrite Eq. (pH)) as 

£ TifCj-i"^ = -Pi 1 '" 11 " + (terms nonlinear in 5/) . (62) 
3=0 

Then we multiply Eqs. ([29|. ([301. and (JHTJ) by 7^, sum over r, and substitute Eq. 1(62]). Next, we use Eqs. ([56} and ([60]) 



to replace all irreducible moments pf 1 ' " Mf appearing in the equations by the fluid-dynamical variables. Additionally, 
all covariant time derivatives of ao, Po, and w M are replaced by spatial gradients of fluid-dynamical variables using 
the conservation laws in the form shown in Eqs. ([55)) . ([M|) . and ([55"|) . The resulting equations of motion are formally 
given as 

mil + n = -ce + j + k + n , 

r„n<"> + n^ 1 = k^Z" + J» + W + W , 
t^ v) + tt^ = Tiqa^ + J^ V + + TO" . (63) 
We remark that in order to derive these equations of motion, it is necessary to use Eq. (|46[) in the following form, 

m 

j=0 Xm 

In the above equations of motion all nonlinear terms and couplings to other currents were collected in the tensors J , 
K, K, W, J^, W, and W" v . The tensors J, J*, and contain all terms of first order in Knudsen and 

inverse Reynolds numbers, 

J = -4inV • n - T nn n ■ F - Snn^O - X nn n ■ I + Xnir^cr^ , 
J» = -n v w vli ~ 5 nn n»9 - 4nV^n + ^ATatt* + TnuTlF" - T n ^ v F v 

- XnuflvCT^ + AnnllJ^ - Xmr^Iv , 

+ 4nV ( "^ > + A„>r>. (65) 

where we defined F^ = W^Pq. In principle, one could replace this quantity by the acceleration u M using Eq. (|35l) . 
The tensors IC, IC^, and KJ 1 " contain all terms of second order in Knudsen number, 



K = Ci + C2 o^ v + Cs 2 + C 4 1 ■ I + Cs F ■ F + Ce I ■ F + (7 V • / + (s V • F, 



+ r/5/ ( "/' y> +f lQ F^F' j) + m I^F v) +r l9 V^F» ) . (66) 
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It is important to remark that among the terms of C(Kn 2 ) is a term uj x ' p oj u "> x . Such a term was believed not to exist 
in a derivation of fluid dynamics from the Boltzmann equation and was therefore speculated to be of quantum nature 
[23| . From our derivation of fluid dynamics, one can see that this is not the case: it simply emerges from a proper 
truncation of the single-particle distribution function. The tensors 1Z, IZ^, and VJ LV contain all terms of second order 
in inverse Reynolds number, 

1Z= tpiU 2 + tp 2 n ■ n + (^371-^71-^, 
TZ" = ip 4 n v TT liV + ^Iln'*, 

1Z^ = ^ 6 nV ll/ + p 7 TT X ^ +ipm {tl n"\ (67) 

In Eq. (62]), terms of order C(Kn 3 ), OiR^RJ 1 ^ 1 ), OiKnR^ 1 ) and 0(KnRr 1 RT 1 ) were omitted. 

Note that we have obtained equations of motion which are closed in terms of 14 dynamical variables. We remark 
that this was accomplished without making use of the 14-moment approximation. This means that the reduction 
of degrees of freedom was not obtained by a direct truncation of the moment expansion, but by a separation of the 
microscopic time scales and the power-counting scheme itself. The information about all other moments are actually 
included in the transport coefficients, as will be shown later. If we also neglect the terms of second order in inverse 
Reynolds number we recover the equations of motion that are of the same form as those derived via the 14-moment 
approximation [20]. However, even in this case, the coefficients in Eqs. (1651) and relaxation times are not the same as 
those calculated from the 14-moment approximation of Israel and Stewart. 

The resulting equations of motion (|63p contain a large number of transport coefficients. In particular, the viscosity 
coefficients and relaxation times of the dissipative currents were found to be, 



n 2 

^ = E^ 2) « r 2) - (68) 

r=0 

Note that in general these transport coefficients depend not only on one moment of the distribution function but on 
all moments of corresponding rank £. As in Chapman-Enskog theory, the viscosity coefficients can only be obtained 
by inverting A^. However, to obtain the transient dynamics of the fluid, characterized by the relaxation times, it is 
also necessary to find the eigenvalues and eigenvectors of AS' ■ 

In practice, the expansion (|18|) is always truncated at some point and the matrices A^ l \ OS 1 ', and t™' will actually 
be finite. The truncation of this expansion was already introduced as an upper limit, Ni, in the corresponding 
summations. In principle, one should only truncate the expansion (|18|) when the values of all relevant transport 
coefficients have converged. Note that different transport coefficients may require a different number of moments to 
converge. 

VII. APPLICATIONS 

In this section, we compute the transport coefficients for several cases. First, we considered the lowest possible 
truncation scheme for Eq. (fl"8]) with Nq = 2, N\ = 1, and N2 — 0. In this case, the distribution function is expanded in 
terms of 14 moments and is actually equivalent to the one obtained via Israel-Stewart's 14-moment ansatz. Second, we 
consider the next simplest case and take No = 3, N% = 2, and -/V2 = 1. Then, the distribution function is characterized 
by 23 moments, and consequently we shall refer to this case as 23-moment approximation. Finally, we include 32 and 
41 moments and verify that the numerical values for the transport coefficients converge. 

We also compute the transport coefficients of the terms appearing in J, J^, and J^ v which are displayed in 
Appendix [C] These transport coefficients were also calculated in previous derivations of fluid dynamics from the 
Boltzmann equation. We shall explicitly point out the corrections to the previous results introduced by our novel 
approach. Note, however, that we are using a linear approximation to the collision term. Nonlinear contributions 
could in principle also enter the transport coefficients in the equations of motion (I63|) , but will not be calculated here. 
Such an investigation will be left for future work. For this reason we also do not compute any coefficient of the terms 
of C^R-iR- 1 ), i.e., entering 1Z, and 1Z^ , since all of them originate exclusively from nonlinear contributions to 
the collision term. 



1 




"lor 




Xo 








T 


E 



(1) ' 

Xo 



r (0) a (0) 

'Or "r i 



(2) ' 

Xo 



r=0,/l,2 



V T (1) aW 

/ j 'Or u r 1 
r=0,#l 
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14-moment approximation 



The 14-moment approximation is recovered by truncating Eq. (|18| at No = 2, N% = 1, and N 2 — 0. For this specific 
truncation A>® is nothing but a number (because for A^ we have to exclude the second and third rows and columns 
and for A^ the second row and column), and thus 



M = 



AW : 



Then, the equations of motion and transport coefficients reduce to those derived in Ref. [2; 

For a classical gas of hard spheres with total cross section <x, in the massless limit, the integrals A^ 1 ' = A, 
A^ = Aq Q can be computed and have the simple form 

A^ = 

A^ = 



(i) 
on 



9A m f p 
3 



5A 



mfp 



and 

(69) 
(70) 



where we defined the mean free-path A m f p = 1/ {tiq<j). The details of this calculation are shown in Appendix FBI 
The coefficients in the ultra-relativistic limit, m/?o — > 0, can then be calculated analytically. The coefficients of order 
©(KnR" 1 ) are collected for the shear stress and particle diffusion in Tables HI and ITll Note that, in this limit, the bulk 



viscous pressure vanishes, and thus we do not need to compute A, 



(0) 





Tn [A mfp] 










Tmr [To,] 


3/ (16a) 


9/4 


1 


3/5 


Pa/20 


Po/20 






TABLE I: The coefficients for the particle diffusion for a classical gas with constant cross section in the ultrarelativistic limit, 
in the 14-moment approximation. 



V 


T"ir[A m f p ] 












4/(3a/? ) 


5/3 


10/7 





4/3 









TABLE II: The coefficients for the shear stress for a classical gas with constant cross section in the ultrarelativistic limit, in 
the 14-moment approximation. 



B. Next correction: 23-moment approximation and beyond 



In order to better understand our formulas, Eqs. (|68[) . we would like to compute the first correction to the expressions 
in Tables U and |ITJ For this purpose, we consider N = 3, Ni = 2, and N 2 = 1. Then, A {1) , fl®, and t& are, after 
removing trivial rows and colums, 2x2 matrices that can be computed from the collision integral Eq. (|45[) . We obtain 
the elements of A^' 2 \ its inverse t^- 1,2 \ and Q^ 1 ' 2 ) as 

AW 1 ( 2 /Sg/30 \ A(2) = 1 ( 9/10 -A/20 \ 



3A m f p 

V "4/3(7 2 1 J ' A mf P I 4 / ( 3 A)) 1/3 

(1) = 3 / 15/4 -ft/s) T (2) = i A f 10 3/V2\ f72) 

1 15/3,, 2 15/2 )' \\ kmiv \-m^ x 27 J' [ ] 

= f / 2 / ^ 2 V fi(2) = ( 1 , / 1 \ (73) 

! (15 + VW5) /3 - 2 (-15 + 7T05) /3 - 2 i I S/3,7 1 lO/S/^ 1 ) 
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see Appendix [Bl for details. The eigenvectors of A^ and A^ are 



X^ 



2A 



(2) 

Xo - 



mfp 

l 



xi 1} 



2\ 



mfp 



(2) 
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(74) 

,x - - r\ ' ^ 

^A m f p tOA m fp 

Using the formulas derived in this paper, Eqs. (|68l) . we calculate the corrected values for the particle- number 
diffusion coefficient and diffusion-relaxation time and for the shear viscosity and shear-relaxation time, 



21 
128 



"oAmfp ^0.164 n A m f p , 



90 



45 - vT05 
14 



Amfp — 2.5897 A m f p 



V = YY P oAmf P ^ 1-2727 P A mfp , 



"7V 2 Amfp , 

where we used that, in the massless and classical limits, 



(76) 
(77) 

(78) 
(79) 



J 1 ) 



„( 2 ) 



1 

12 



n , a 2 



(i) 



Ad 



Po, 



„( 2 ) 



0o 



Po- 



As before, the coefficients in the ultra-relativistic limit, m/3o — > 0, can then be calculated analytically. The coefficients 
of order ©(KnR" 1 ) are collected for the shear stress and particle diffusion in Tables UTT1 and HVl 









Ann [Tto] 


An7r \y~ri\ 




"7~n TV \p~n ] 


21/(128o-) 


2.59 


1.00 


0.96 


0.054/3 


0.118/30 


0.0295/3 /P 



TABLE III: The coefficients for the particle diffusion for a classical gas with constant cross section in the ultrarelativistic limit, 
in the 23-moment approximation. 



V 


1"7r[A m f p ] 




A 7T n ["TV ] 








14/(11(7/30) 


2 


134/77 


0.344/3(7 1 


4/3 


-0.689/3(7 1 


-0.689/no 



TABLE IV: The coefficients for the shear stress for a classical gas with constant cross section in the ultrarelativistic limit, in 
the 23-moment approximation. 

To obtain these expressions we used the results from Appendix [Dl and that, in the massless/classical limits, D20 = 
3Pq . Note that most of the transport coefficients were corrected by the inclusion of more moments in the computation. 
The coefficients related to the shear-stress tensor were less affected by the additional moments, when compared to the 
particle-diffusion coefficients. This might explain the poor agreement between the Israel-Stewart theory and numerical 
solutions of the Boltzmann equation in Refs. [17j regarding heat flow and fugacity. 

We further checked the convergence of this approach by taking 32 and 41 moments. In this case, the matrices A^ 1,2 \ 
t' 1,2 ) and f^ 1 ' 2 ) were computed numerically. There is a clear tendency of convergence as we increase the number of 
moments. For the particular case of classical particles with constant cross sections, 32 moments seems sufficient. See 
Tables |V] and Ell for the results. 



VIII. DISCUSSION AND CONCLUSIONS 



A. Knudsen number and the reduction of dynamical variables 

It is important to mention that the terms IC, KJ 1 , and K,^ which are of second order in Knudsen number lead to 
several problems. The terms which contain second-order spatial derivatives of u M , ao, and Pq, e.g., V M / AI , V^F^, 
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number of moments 




Tn [A m fp] 




Ann \p~n\ 


An7T [7"n] 






14 


3/(16<x) 


9/4 


l 


3/5 


ft/20 


ft/20 





23 


21/ (128cr) 


2.59 


1.0 


0.96 


0.054ft 


0.118ft 


0.0295ft /Po 


32 


0.1605/a 


2.57 


1.0 


0.93 


0.052ft 


0.119ft 


0.0297ft /Po 


41 


0.1596/a 


2.57 


1.0 


0.92 


0.052ft 


0.119ft 


0.0297ft /Po 



TABLE V: The coefficients for the particle diffusion for a classical gas with constant cross section in the ultrarelativistic limit, 
in the 14, 23, 32 and 41-moment approximation. 



number of moments 


n 


Ttt [A m f p ] 




A TYn [Ttj- ] 








14 


4/(3<rft) 


5/3 


10/7 





4/3 








23 


14/(llaft) 


2 


134/77 


0.344ft 1 


4/3 


-0.689/ft 


-0.689/no 


32 


1.268/(crft) 


2 


1.69 


0.254ft" 1 


4/3 


-0.687/ft 


-0.687/no 


41 


1.267/(crft) 


2 


1.69 


0.244ft 1 


4/3 


-0.685/ft 


-0.685/no 



TABLE VI: The coefficients for the shear stress for a classical gas with constant cross section in the ultrarelativistic limit, in 
the 14, 23, 32 and 41-moment approximation. 

V<"/">, \7^F"\ A£d„<T a ", and Vfl, are especially problematic since they change the boundary conditions of the 
equations. In relativistic systems these derivatives, even though they are space-like, also contain time derivatives and 
thus require initial values. This means that, by including them, one would have to specify not only the initial spatial 
distribution of the fluid-dynamical variables but also the spatial distribution of their time derivatives. In practice, 
this implies that we would be increasing the number of fluid-dynamical degrees of freedom. 

There is an even more serious problem. By including terms of order higher than one in Knudsen number, the 
transport equations become parabolic. In a relativistic theory, this comes with disastrous consequences since the 
solutions are acausal and consequently unstable Q . For this reason, if one wants to include terms of higher order in 
Knudsen number, it is mandatory to include also second-order co-moving time derivatives of the dissipative quantities. 
Or, equivalently, one could promote the moments p^tPi or further ones to dynamical variables. For this reason 
we do not compute the transport coefficients for these higher-order terms in this paper. 

In practice, a way around this would be to replace e.g. the a x ^a^ term in K^ v using the asymptotic (Navier-Stokes) 

solution by (I^tt^ct^, and thus effectively rendering it a term contributing to J^ v . This should be a reasonable 
approximation if one is sufficiently close to the asymptotic solution. This would then change the coefficient of the 
respective term in J^ v . In principle, this could be done to all terms in IC, /C M , and TO 1 " , except for the ones containing 
exclusively powers and/or gradients of and w^. In the same spirit, using the asymptotic solutions one could also 
shuffle some of the terms in J , , and J^ u (those not containing F^, w* 1 ", and gradients of dissipative currents) 
into terms contributing to 1Z, TV 1 , and TV 1 " (or vice versa). How this changes the actual transient dynamics remains 
to be investigated in the future. 

B. Navier-Stokes limit 

Note that one of the main features of transient theories of fluid dynamics is the relaxation of the dissipative currents 
towards their Navier-Stokes values, on time scales given by the transport coefficients Tn, t„, and r n . From the 
Boltzmann equation, Navier-Stokes theory is obtained by means of the Chapman-Enskog expansion which describes 
an asymptotic solution of the single-particle distribution function. It is already clear from the previous section that 
the equations of motion derived in this paper approach Navier-Stokes-type solutions at asymptotically long times, in 
which the dissipative currents are solely expressed in terms of gradients of fluid-dynamical variables. 

It is interesting to investigate, however, if our equations approach the correct Navier-Stokes theory, i.e., if the 
viscosity coefficients obtained via our method are equivalent to the ones obtained via Chapman-Enskog theory. It 
should be noted that this is not the case for Grad's and Israel and Stewart's theories [a [HI H3- The viscosity 
coefficients computed by these theories do not coincide with those extracted from the Chapman-Enskog theory. We 
remark that, after taking into account the first corrections to the shear viscosity coefficient, see Eq. (l78|) and Table IVTl 
our result approached the solution obtained using Chapman-Enskog theory, ijns = 1-2654/ (/3o<j) @- In principle there 
is no reason for the method of moments to attain a different Navier-Stokes limit than Chapman-Enskog theory. We 
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can show that, if the same basis of irreducible tensors k^ 1 ■ ■ ■ and polynomials P\£ is used in both calculations, 
they both yield the same result, even order by order. 

C. "Non-hydrodynamic" modes and the microscopic origin of the relaxation time 

One of the features of the theory derived in this paper (and also of Grad's and Israel-Stewart's theories) is the 
appearance of so-called non- hydrodynamics modes, i.e., modes that do not vanish in the limit of zero wave- number. 
Such modes do not exist in Navier-Stokes theory or its extensions via the Chapman-Enskog expansion. For this 
reason, these modes are usually not associated with fluid-dynamical behavior, hence the label "non-hydrodynamic". 

The non-hydrodynamic modes describe the relaxation of the dissipative currents towards their respective Navier- 
Stokes solutions and can be directly related to the respective relaxation times. For the case of the shear non- 
hydrodynamic mode, w"[^" a ~ hydro (k), it can be shown that in the limit of k — > the mode is given by w"™ a ~ hydl ° (0) = 
— i/Tn [5j]. In Chapman-Enskog theory the transient dynamics of the system is neglected, e.g., it is assumed that in 
the absence of space-like gradients, time-like gradients vanish as well, and it is natural that such modes do not exist. 

The appearance of non-hydrodynamic modes in a fluid-dynamical theory seems to counteract the prevalent belief 
that fluid dynamics effectively describes the asymptotic long-time and long-distance behavior of the microscopic 
theory. Recently, a microscopic formula for the relaxation time of dissipative currents was obtained in the framework 
of linear response theory |10| . In that paper, the relaxation time was shown to be intrinsically related to the slowest 
microscopic time scale of the system, i.e., to the singularity of the retarded Green's function closest to the origin in 
the complex-plane. Thus, the non-hydrodynamic modes in Israel and Stewart's theory and in the equations derived 
in this paper belong to a description at long, but not asymptotically long, times. 

This means that the theory derived in this paper (as well as Israel and Stewart's theory) attempts to describe the 
dynamics of the dissipative currents at time scales of the order of the (slowest) microscopic times scale (which is of 
the order of the mean-free path). Such findings challenge the point of view that a fluid-dynamical description can 
only be formulated around zero frequency and wave number and that the inclusion of relaxation times can only be 
understood as a regularization method to control the instabilities of the gradient expansion. In fact, the relaxation 
times correspond to microscopic time scales, independent of any macroscopic scale related to the gradients of fluid- 
dynamical variables. Note that the expressions presented in Ref. [To| and in this paper for r\ and r w are equivalent. 

D. Conclusions 

In this work we have presented a general and consistent derivation of relativistic fluid dynamics from the Boltzmann 
equation using the method of moments. First, a general expansion of the single-particle distribution function in terms 
of its moments was introduced in Sec. IIIII We constructed an orthonormal basis which allowed us to expand and 
obtain exact relations between the expansion parameters and irreducible moments of the deviations of the distribution 
function from equilibrium. We then proceeded to derive exact equations for these moments. 

The main difference of our approach to previous work is that we did not close the fluid-dynamical equations of 
motion by truncating the expansion of the distribution function. Instead, we kept all terms in the moment expansion 
and truncated the exact equations of motion according to a power-counting scheme in Knudsen and inverse Reynolds 
number. Contrary to many calculations, we did not assume that the inverse Reynolds and Knudsen numbers are of 
the same order. As a matter of fact, in order to obtain relaxation-type equations, we had to explicitly include the 
slowest microscopic time scales, which are shown to be the characteristic times within which dissipative currents relax 
towards their asymptotic Navier-Stokes solutions. Thus, Navier-Stokes theory, or the Chapman-Enskog expansion, is 
already included in our formulation as an asymptotic limit of the dynamical equations derived in this paper. 

We concluded that the equations of motion can be closed in terms of only 14 dynamical variables, as long as we only 
keep terms of second order in Knudsen and/or inverse Reynolds number. Even though the equations of motion are 
closed in terms of these 14 fields, the transport coefficients carry information about all moments of the distribution 
function (all the different relaxation scales of the irreducible moments). The bulk- viscosity, particle-diffusion, and 
shear-viscosity coefficients agree with the values obtained via Chapman-Enskog theory. 
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Appendix A: Derivation of the collision terms 

In this appendix, we derive Eqs. (pH)) and P5|) . The first step is to linearize the collision operator, 

C{f] = lJ dK'dPdP'W^ pp , (/p/p'/ k / k < - /k/k'/p/p') , (Al) 

in the deviations from the equilibrium distribution functions. In the main text, the deviations from the local- 
equilibrium distribution function were parametrized as 

$fp = f P - fo P = /op/op^p • ( A2 ) 
Then, only keeping terms of first order in <j>, we can prove that 

/p/p< - fop fop- (l + fo P 4p> + /o P 0p) + O {4?) , (A3) 
/p/p' = /op/op'(l-a/o P '0p'-a/op</' P )+O(0 2 ) . (A4) 
Substituting Eqs. (|A"3t and (|A4| into Eq. (TAT]) , we obtain, 

C[f] = ^f dK'dPdP'W^ pp J ufovfophp> {4>p + 4>p' + {4> 2 ) , (A5) 

where we also used the equalities 

fop = fo P exp((3 E p - a ) , (A6) 

/op/op'/ok/ok' = /ok/ok'/op/op' • (A7) 

Inserting Eq. (|A5|) in the expression for the irreducible collision term (f3l?|) , we obtain 

C^Y" W> = \ J dKdK'dPdP'W^ ppl f ok f vf J 0p , 

x^ 1 ^" 1 • • • fcW> (^p + <t>V - fa - <h>) + O (0 2 ) . (A8) 

The next step is to substitute the moment expansion of the single-particle distribution function, Eqs. (fTTj) and (|24p . 
into Eq. (|A8[) . expressing it in the following form 

oo N m 

C^V W> = - E E (ArrXPZ #""*" + (f) > ( A9 ) 

7/7=0 71— 

where we defined the tensor 

(A n ):;::* = \ j dKdK'dPdP'W^^rfoufovfoJop'EZr 1 ^ •••fc"<> 

x (n™k {vi ■■■k„ m) +H { ^k' {l/i ...k' Vm) -n^p {Vl ... Pum) -n^ P ' (ui ... P ' Um) ) . (aio) 

The integral (An)^ 1 ...^' is a tensor of rank m + £, which is symmetric under permutations of /z-type indices and 
symmetric under permutations of v—type indices, and which depends only on equilibrium distribution functions. The 
latter contain only the fluid 4-velocity u^. Therefore, (Am)^...^" must be constructed from tensor structures made 
of and the metric tensor g^ v . Also, (An)^ 1 .'.'.'^ was constructed to be orthogonal to u M and to satisfy the following 
property, 

a «'"W A-An \-^rn) Ul . .. Vm ~ K^rn) f 3l ...p m ■ 

Since {A rn )^l^ e is orthogonal to u M , it can only be constructed from combinations of projection operators, . 
This already constrains m + £ to be an even number, since it is impossible to construct odd-ranked tensors solely from 
A^s. This means that both £ and m are either even or odd. Therefore, the following type of terms could appear in 

I A M/ . 

\S^rn) v ...„ • 
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(i) Terms where all /i— type indices pair up on projectors A Mi '^ and all ^-type indices on projectors A„ „ , e.g. 

A ^M2 . . . AWW . . . A^-^'A^ • • • A„ p „ a • • • A Vm _ lUm . (A12) 

All possible permutations of the //-type indices among themselves and ^-type indices among themselves are 
allowed. 

(ii) Terms where at least one //-type index pairs with a f-type index on a projector, e.g. 

A£A" a " 3 ...A^ ■■•A w - lW A^ •••A VpV , •••A„ m _ 1 „ m . (A13) 

Again, all possible permutations of the //-type and i/-type indices are allowed. If there is an odd number of 
projectors of the type A£*, both £ and m must be odd. If there is an even number, both £ and m must be even, 

too. Without loss of generality, suppose that t > m. For I + m to be even, I must be m + 2, m + 4, Then 

one could pair all i/— type indices with //-type indices on projectors of the form A^*, with some projectors left 
over which carry only //-type indices, e.g. AW*. 

(iii) If ^ = m, all //-type indices could be paired up with z/-type indices on projectors of the form A^ , with no 
left-over projectors like explained at the end of (ii), 

A£..-A«. (A14) 

Again, all permutations of the //-type indices among themselves and v—type indices among themselves are 
allowed. 

Note that terms of the type (i) and (ii) by themselves do not satisfy the property (|A11[) . This happens because any 
term which contains at least one projector of the type A MiAtj or A VpVq vanishes when contracted with A^.'.'.'^ A^...^™ . 
Therefore, (An)^ 1 .'.'.^ cannot be solely constructed from terms of type (i) and (ii), because otherwise it would vanish 
trivially, and property (IA11|) would not be satisfied. There must at least be one term of type (iii). However, this 
implies that m = £. This does not imply that terms of type (i) and (ii) do not appear; they do occur, but in such a 
way that Eq. (|All|) is satisfied. In summary, {Am)^...^ has the form 

(A„CX = 6 tm {A%A$ ■ • • A + [terms of type (i) and (ii)] } , (A15) 

where the parentheses denote the symmetrization of all Lorentz indices. Contracting Eq. (IA15|) with A^J.'.'.'^A^ " 
and using Eq. (|Alip , we prove that 

{Am)t;::z = <w 41 3 ^f'Z • (Aie) 

Finally, substituting Eq. (IA16|) into Eq. (|A9|) we derive Eq. (|3I|), introduced in the main text of the paper, 

oo 

q£r w> = -£^rfr w - ( Ai7 ) 

m— 

The coefficients Arn can be obtained from the following projection of (Am) 1 ^ 1 ' 



J\W — A" 1 "'"' (A 

— fit 



= - ( J + t) y dKdK'dPdP'W kk ^ pp ,f a U bk '/op i "op' Bt -1 *^ 1 •••fc w> 

X (^ki *<Mi ' ' ' fc M + ^k'l ^ ' ' ' Ke) - P<Mi <l P<Mi ' ' ' * «>) > ( Al8 ) 

where we used that A^"!{^ = 2£ + 1. 
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Appendix B: Calculation of the collision integrals 



In this appendix, we calculate the collision integrals, Eq. (|45|) . for a classical gas, i.e., fa k = 1, of hard spheres in 
the ultrarelativistic limit, m/3g <§C 1. Then, Eq. (|45l) becomes 



A {e) 



v{2£ 



l—j J dKdK'dPdP'W kkl ^ pp ,f ok f okl 



x E^k^ ...*•*> (nilk (ui ■ ■ ■ k v>) + n$ n k' {ui ■ ■ ■ K t) - n p e lp Wl ■ ■ -p Vl) - n% P ' {ui ■ ■ -p' Ve) ) . (bi) 

The functions n k ^ were defined in the main text, see Eq. (|25[) . The transition rate Wkk/->pp/ is written in terms of 
the differential cross section a(s, 0) as 

WWpp, - sa(s, 6) (2ir) 6 6^ (k^ + fc'" - pf* - p'») . (B2) 
The variable s and are defined as 

S = {k + k')\ C osQ= {k - k ' ) - (p : p,) . (B3) 

{k-k'f 

We further define the total cross section as the integral 

a T (s) = ^ J d<d sin9CT(s,e) . (B4) 

In order to calculate Am it is convenient to first define the tensors X™. „ 

(lis yi / m 

X£ v „ m = i J dKdK'dPdP'W kk ,^ ppf f ok f k' 

x E^k^k,; (fc 71 • • • fc 7m + fc 7i • • • k!y m — Pj 1 ■ ■ ■ Pj m — p 7l • • ' p'-y m ) ; (B5) 

The collision integrals Am can always be expressed as linear combinations of contractions/projections of -^uw 7l ... 7m - 
For the purpose of this paper, we shall only need Xuwy 1 —-y m for m = 2 and 3. For now we concentrate on calculating 
these integrals. We separate X™ as 

yn _ An , on fRfi") 

^^7l"'77Ti A i ^7i-"7?ji M l/ 7l"'7m ' V / 

with 

^ 7l ... 7m = ~ / dKdK'dPdP'W kk ^ pp ,f ok fowE^ k^K (fc 71 • • • fc 7m + fc 7l • • • k' lm ) , 

S^ 7l ... 7m = ~ j dKdK'dPdP'W kkl ^ ppl f ok f ok >EZ k u k„ (p 7l • • -p Jm + p 71 • • • p' lm ) . (B7) 

The dPdP' integration in the first tensor, A™ , can be immediately performed and written in terms of the total 

cross section, ctt(s), as 

^ 7l ... 7m = / dKdK'f ok fo k 'EZk u k„ (fc 7l • • • k lm + k' lx ■ ■ ■ k' lm ) S ~G T (s) . (B8) 
The calculation of the second tensor, B" 7l „ , is cumbersome. First, we write it in the general form 



TYn 



j dKdK' '/ok/ok'-Ep^e-n...-^, (B9) 



where we introduced the tensor 



© 7l - 7m = \ J dPdP'W kk ,^ pp , Pll • • -p 7m . (B10) 

The integral © 7l ... 7m is an m-th rank tensor. Strictly speaking, for isotropic cross sections, this tensor can only depend 
on the normalized total momentum of the collision = s" 1 ^ 2 (k^ + k'^) = s _1 / 2 P^ . Thus, the tensor structure of 
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7l ... 7m must be constructed by combinations of P^ and the projection operator orthogonal to P^, Ap" = g^ v — P^PJf. 
In general, 

[m/2] 

®7l— 7m = ^ ' ("I)* a mqCmqCS /1 ... 7m , (Bll) 
9=0 

where we defined 

(2g-l)H, 



(m - 2q)\2q\ 



CI _ A (7172 A 72,-172, p72,+ l 57m) 

Cm, = v{2q+l y A J dPdP'W kk ^ pp , [P^j (-Afpapp) . (B12) 

The parentheses () denote the symmetrization of the tensor. For example, 

6 7l72 = C20-PT71-PT72 — C2iAp 7l72 , 
©717273 — C3qPt 1i Pt 12 Pt 13 — C31 ^Ap 7i72 Pt 73 + A Plll3 P Tl2 + Ap 7273 P T 7i^ , (B13) 

The integrals C nq are scalars and can be computed in any frame. It is most convenient to calculate them in the 
center-of-momentum frame, where, Pj, = (1,0,0,0) and = diag (0, — 1, — 1, — 1). Then, it is straightforward to 
prove that 

C nq = f T (s) ein-zo+W ( s - 4m 2 ) (2g+1)/2 = ^(s) {n+2)/2 ( } 

In the massless limit, the tensors Jf™ and X™„ 7i7273 become 

^7i7 2 = J dKdK' fokfok'K^Ko-T (s) k x k' x (k 7l k 12 + k'^k'^ - ^P Tii Pt J2 + ^s 57l72 ^ , 
^7!7273 = / dKdK'fo^E^Kar (s) k x k' x 

7172 Pt-Is + 07173-^72 + 37273 Prii) 



fc 71 k 12 k~ (3 + fc 7i fc 72 k l3 — Pru Pt~!2 Pt^ 3 



(B15) 



where we used that, in the massless limit, s = 2k x k' x . 



1. Particle-diffusion current 

For the collision integrals related to the particle-number diffusion current, we need the following two contractions 

A^U v U^X^ 7ll2 = -U T (hoIn+5,1 - 4/2l/„+4,l - / 3 l/n+3,l) , 
A«%V^i; 7l7J73 = (37 10 /n+6,l - H/2l/„+5,l " 5/ 3 i7 n+4) l - 3/ 4 lin+3,l) ■ (B16) 



To obtain the above relations, we used Eq. (fT6f and the definitions (]39|). In the massless and classical limits the 
integrals I nq = J nq can be calculated analytically 

" 9 ■ 9 (2 g +l)!!2^ /3 "+ 2 (2<7+l)!!2^- 2 ' 1 ' 
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Then, 



A* W71 72 



4 P 



A'< 7 W 2 X°,, 7l72 =-24n a T §, 

Po 



A M7l w' y U 72 W 73 X~ i = 12?ln(7T — - 



pT 

/3o 4 ' 



A" 7l ^ u 72 M 73 A^ 7i7273 = -280n ot T T -^. (B18) 



As a consistency check, we confirmed that A'' 7i m 1/ !i 72 X (1 J 7i72 = A /i7l u I/ u 72 u 73 X /xl/ 1 7i7273 = 0. 

The components of A^' change according to the number of moments included. In the 14-moment approximation, 
using Eqs. (|2"5|) and (flU)) . we obtain 

4o> = ^afoMl'A^/^I^^ = \n a T . (B19) 
In the 23-moment approximation, e.g. considering three polynomials in the expansion (fT8|) . for I = 1, 

.4 



(1) 


_ W« 


rO 


3 


(1) 




r2 


3 



71 72 73 

= — ( »^«-;; , ^ ; 1 « v ^. v:,: ; i .^'4;'a«/^^i;-; W3 j . (B20) 

Then, using the results from Appendix [El for the coefficients ai\ together with Eqs. (|B17I) and (|B18|) . we obtain 

/i(i) 2 /i(i) Po 

(1) ^ fl) 1 

^20 = ~i^2 n 0O"T, ^22 = q n oor- (B21) 



2. Shear-stress tensor 



For the collision integrals related to the shear-stress tensor, we need the following two contractions 

A^ 7l72 A^ 7i72 = ya T (hoIn+5,2 + 4/ 2 l/„ + 4, 2 ) , 
A^ 7l72 W 73 A;„ 7i7273 = 5a T (7 10 I„ +6 ,2 - Wn+5,2 + 2I 31 I„+4,2) ■ (B22) 
In order to obtain the above relations, we used Eq. ([T6|) and the definitions (|39|) . Using Eq. ()B17|) . 
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A^u^X-^ 7i = 132a T J, 

A^u-ri^ = 880a T ^. (B23) 

Po 

The components of A^ change according to the number of moments included. In the 14-moment approximation, 
using Eqs. ((25)) and (p~9|) . we obtain 



4? = ^?T A " l ' 7l72 ^7 1 72 = 5«oa T , (B24) 



22 



where we used Eqs. (IB17|) and (|B23|) . together with the results from Appendix [El 

In the 23-moment approximation, e.g. considering two polynomials in the expansion (|18[) . for £ = 2, 



/|( 2 ) — _ ( 1 i n ^„^\ A A^7l72 V-r-1 , „( 2 )„( 2 ) A/f 7i72„73 

~~ 10 V / ^V"7i72 + "10 "11 " " ^-pi'717273' 

-^rl — "11 "10 ^ ^A"^7l72 + JO 11 11 ^^^717273- 

Then, using once more the results from Appendix IE1 and Eqs. (|B17|) and (|B23|) . we obtain 



f 2 ) ^ (2) 

A w = ^ a T , 4? = i^T- (B26) 

We did not calculate the coefficients related to the bulk viscous pressure, since this quantity vanishes in the massless 
limit. Also, if the mass was taken to be finite, some of the steps taken in this appendix would not be possible. 



Appendix C: Transport coefficients 

In this appendix we list all the transport coefficients of fluid dynamics calculated in this paper. The transport 
coefficients for the bulk viscous pressure are 
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where ho — (so + Po)/™o is the enthalpy per particle. The transport coefficients for the particle-diffusion current 
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The transport coefficients for the shear-stress tensor are 
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Appendix D: Calculations 

(2) 

In this appendix we compute the quantity 7} in the 14-moment approximation and the 23-moment approximation. 
This variable was defined in the main text and is given by, 
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The first step is to compute the thermodynamic integral Trn , 
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In this case, N\ = 1 and A 2 = 0, and 



Also, in the 14-moment approximation, 



1. 14-moment approximation 



~(2) _ T-(2) 
71 — -Mo • 



(D3) 



In the massless/classical limits 



and finally 



y(2)_ (2) p (2)_ W( 2 ) 

k0 — 2! 00 — 2! 



oy(2) _ £0 
" k0 " 8iV 



(D4) 



2. 23-moment approximation 



In this case, Ni = 2 and A 2 = 1, and 



(2) _ (2) n (2) ^(2) 
71 — -MO + "10 -Ml • 



(D5) 



Also, in the 23-moment approximation. 



0,(2) 
rt k0 



2! Mo<> = 

^ (2) J2) p (2)_ ^(2) 



2! 



(ag)) a + ^a<K' 



-1/(2) _ ■ < pi - 1 _ 

kl ~~ 2! 11 1 — 



2! 



(D6) 
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We know that 



Thus, 



and 



(2) 



4Po' y ail > ~ 6' a f)~ ft ■ (D7) 



W £ ) = ^(- 1 + ^k)- (DM 



^ = S^fdK /„ k £ k 1 (-1 + W) (A^Kkp) 2 = -|. (D9) 



4P 5!! y JUK k V 6™ V v a py 30 
Substituting fi( 2 ) from Eq. ([711) we obtain 



7p ) = ^/?o = 0.133A). (D10) 



Appendix E: Orthogonal Polynomials 

In this appendix, we construct the set of orthogonal polynomials used in the main text. These will be polynomials 
in energy, E^ = u M fc M , i.e., orthogonal polynomials generated by the set 1, E^, . . .. We construct this orthogonal 
set using the Gram-Schmidt orthogonalization method. First we introduce 



W = 



(2 , +T yjT (A Q/3 fcaM /ok/ok, (El) 

where / k is the equilibrium distribution function as defined in the main text. The weight W^' will be determined 
such that the orthogonal polynomial of order n = and index £ is normalized, 

dK^P^pW = l, (E2) 

Without loss of generality, the polynomials of order are set to 1 for all values of £, 

P$ = = 1 . (E3) 
Then the normalization parameter is obtained from Eq. (|E2|) . 

W w = {~l) e — . (E4) 
■ht.i 

The thermodynamic functions J nq were defined in the main text, see Eq. (1391) . 
The polynomials are parametrized as 

n 



We construct the polynomials in sequence according to the parametrization (|E5[) starting from n — 0, Eq. (|E3|) . using 
the orthonormality condition (|20l) . The orthogonality /normalization condition implies that, for a polynomial of order 



i P (i) 



20 



for all j < i. Substituting Eq. (|E5|) . we obtain the following equation for the coefficients 

a (i) T 

V V iU) 13 - ■ X, - (V7) 
l^ U kj -77) - , „^2 d >»> l^O 



where k = 0, . . . ,i, and we defined the (i + 1) x (i + 1) matrix T>^ = Jk+j+2£,£- The solution of Eq. (|E7I) is 

(«<?) 2 = (s- 1 )<"^«. 

where (l? -1 )^*'' is the inverse of D^ 1 '. For example, for any polynomial of order 1, the coefficients are 

J2£+2,eJ2£,£ — (J2£+l,l) 2 

(E9) 



l 10 _ \ U J 01 _ J 2£+l 



Appendix F: Irreducible tensors 

In this appendix, we give some practical relations concerning the irreducible tensors k^ 1 k^ 2 ■ ■ ■ fe w ' introduced in 
the main text. The definition of these tensors is 

fc <Pi fc M2 . . . fc w> = ^w-w^i ...h vi . (Fl) 
The projection operator A^j.'.'.'/f/ is symmetric and traceless in the indexes /i and v 

Wi A£;;:# = A£;;;« =0, V l < i, i < I (F2) 

The parentheses () on the indices denotes symmetrization of the tensor. These projections are constructed in Ref. @ 
and can be obtained from 

[1/2] 
fc=0 

CUk)- (-l) k {£l)2 {2£ - 2k)l (F3) 



where in the last summation the symbol [£/2] denotes the largest integer not exceeding t/2 and 

/ 2 fe \ 2 

®(iie) = _ 2k)\ [ — — - J ^Mi/*a . . . ^A»2fe-i/i2fc . . . £Vzk-iV2k ^Atak+i^fc+i . . . (F4) 



This summation is supposed to run over all distinct permutations of ^t-type and i/— type indices (we do not permute 
the indices /i with z/). For I — 2 this recipe gives the usual double symmetric and traceless projection operator 
A^g commonly employed in relativistic fluid dynamics. As mentioned in the main text, this set of tensors are useful 
because they form an orthogonal basis, see Eq. (fTtj)) . 
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